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Abstract 


In  this  report  we  document  the  implementation  of  high-order  Higdon  nonreflecting 
boundary  conditions.  We  suggest  a  way  to  choose  the  parameters  and  demonstrate 
numerically  the  efficiency  of  our  choice.  The  model  we  used  is  the  linearized  2-D  Euler 
equations  with  zero  advection.  These  equations  are  solved  by  the  finite  difference 
method.  We  close  with  a  list  of  topics  for  research. 


1.  Statement  of  the  Problem 

Consider  the  linearized  Euler  equations  in  an  infinite  domain.  For  simplicity  we  assume 
that  the  domain  has  a  flat  bottom  and  that  there  is  no  advection  and  no  Coriolis  effect, 
although  these  assumptions  may  be  removed  in  future  studies.  A  Cartesian  coordinate 
system  (x,  y )  is  introduced,  as  shown  in  the  figure. 

y 


Figure  1:  An  infinite  domain 
The  linearized  Euler  equations  are  (see,  e.g.,  [32]  or  [33]): 


dtp  +  p0(dxu  +  dyv)  =  0  ,  (1) 

dtu+—dxp  =  0,  (2) 

Po 

dtv  +  —dyp  =  0  ,  (3) 

Po 

dtP  +  'yPoidxii  +  dyv)  =  0  .  (4) 


Here  t  is  time,  u(x,  y,  t)  and  v(x,  y,  t)  are  the  unknown  velocities  in  the  x  and  y  direc¬ 
tions,  p(x,y,t)  is  the  density,  p(x,y,t)  is  the  pressure  and  7  is  the  gas  constant.  The 
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linearization  was  done  about  mean  zero  velocities,  mean  density  po  and  mean  pressure 
Po-  We  use  the  following  shorthand  for  partial  derivatives 


dl 


&- 

da1 


The  nonlinear  Euler  equations  are 


dtp  +  dx(pu )  +  dy(pv)  =  0  , 

(5) 

dtu  +  udxu  +  vdvu  -| — dxp  =  fv  , 

P 

(6) 

dtv  +  udxv  +  vdvv  H — dvp  =  —fu  , 

p 

(7) 

dtp  +  udxp  +  vdyp  +  7 p(dxu  +  dyv)  =  0  . 

(8) 

It  can  be  shown  that  a  single  boundary  condition  must  be  imposed  along  the  entire 
boundary  to  obtain  a  well-posed  problem.  At  x  — >  oo  the  solution  is  known  to  be 
bounded  and  not  to  include  any  incoming  waves.  To  complete  the  statement  of  the 
problem,  initial  values  for  u,  v,  p  and  p  are  given  at  time  t  =  0  in  the  entire  domain. 

We  now  truncate  the  infinite  domain  by  introducing  an  artificial  east  boundary  T e, 
located  at  x  =  sg,  and  similarly  at  the  the  other  three  sides  (see  Figure  1).  To  obtain 
a  well-posed  problem  in  the  finite  domain  If  we  need  a  single  boundary  condition  on 
each  of  the  artificial  boundaries  Te,  T\y,  Tat,  and  Tg.  This  should  be  a  Non-Reflecting 
Boundary  Condition  (NRBC).  We  shall  apply  a  high-order  NRBC  for  the  variables.  A 
discussion  on  this  NRBC  follows. 


2.  Higdon’s  NRBCs 

On  the  artificial  boundary  T £  we  use  one  of  the  Higdon  NRBCs  [1].  Similarly  for 
the  other  three  sides.  These  NRBCs  were  presented  and  analyzed  in  a  sequence  of 
papers  [2] — [6]  for  non-dispersive  acoustic  and  elastic  waves,  and  were  extended  in  [1] 
for  dispersive  waves.  Their  main  advantages  are  as  follows: 

1.  The  Higdon  NRBCs  are  very  general,  namely  they  apply  to  a  variety  of  wave 
problems,  in  one,  two,  and  three  dimensions  and  in  various  configurations. 

2.  They  form  a  sequence  of  NRBCs  of  increasing  order.  This  enables  one,  in  principle 
(leaving  implementational  issues  aside  for  the  moment),  to  obtain  solutions  with 
unlimited  accuracy. 

3.  The  Higdon  NRBCs  can  be  used,  without  any  difficulty,  for  dispersive  wave  prob¬ 
lems  and  for  problems  with  layers.  Most  other  available  NRBCs  are  either  de¬ 
signed  for  non-dispersive  media  (as  in  acoustics  and  electromagnetics)  or  are  of 
low  order  (as  in  meteorology  and  oceanography). 

4.  For  certain  choices  of  the  parameters,  the  Higdon  NRBCs  are  equivalent  to  NR¬ 
BCs  that  are  derived  from  rational  approximation  of  the  dispersion  relation  (the 
Engquist-Majda  conditions  being  the  most  well-known  example).  This  has  been 
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proved  by  Higdon  in  [1]  and  in  earlier  papers.  Thus,  the  Higdon  NRBCs  can  be 
viewed  as  generalization  of  rational-approximation  NRBCs. 

The  scheme  developed  here  is  different  than  the  original  Higdon  scheme  [1]  in  the 
following  ways: 

1.  The  discrete  Higdon  conditions  were  developed  in  the  literature  up  to  third  order 
only,  because  of  their  algebraic  complexity  which  increases  rapidly  with  the  order. 
Here  we  show  how  to  easily  implement  these  conditions  to  an  arbitrarily  high 
order.  The  scheme  is  coded  once  and  for  all  for  any  order;  the  order  of  the 
scheme  is  simply  an  input  parameter. 

2.  The  original  Higdon  conditions  were  applied  to  the  Klein-Gordon  linear  wave 
equation  and  to  the  elastic  equations.  Here  we  show  how  to  apply  them  to  the 
linearized  Euler  equations  (l)-(4). 

3.  The  Higdon  NRBCs  involve  some  parameters  which  must  be  chosen.  Higdon  [1] 
discusses  some  general  guidelines  for  their  manual  a-priori  choice  by  the  user.  We 
shall  show  how  a  simple  choice  for  these  parameters  can  dramatically  simplify 
the  calculations  and  enable  implementation  of  NRBCs  of  much  higher  order  with 
less  computational  overhead. 

The  Higdon  NRBC  of  order  J  is 


Hj  : 


n  ($ + 


3= 1 


i)  =  0  on  Te 


(9) 


Here,  the  Cj  are  parameters  which  have  to  be  chosen  and  which  signify  phase  speeds 
in  the  x-direction.  The  boundary  condition  (9)  is  exact  for  all  waves  that  propagate 
with  an  x-direction  phase  speed  equal  to  any  of  C\ ,  . . . ,  Cj.  This  is  easy  to  see  from 
the  reflection  coefficient  (see  below).  For  the  other  sides  we  replace  dx  by  the  normal 
derivative  to  the  boundary. 

We  make  a  few  observations: 

•  In  their  sequence  of  papers  [22] — [31],  Givoli,  Neta,  and  van  Joolen  showed  that 
one  should  always  take  Cj  >  Co  because,  in  general,  the  solution  consists  of  an 
infinite  number  of  waves  with  different  phase  speeds.  For  this  problem,  however, 
there  is  only  one  wave  speed  Co  impinging  on  the  boundary  at  all  possible  angles; 
hence,  one  should  instead  take  Cj  <  Co  (i.e.,  Cj  =  Co  cos (<j>j),  where  (f>j  GO...  |). 

•  The  first-order  condition  H\  is  a  Sommerfeld-like  boundary  condition.  If  we  set 
C\  =  Co  we  get  the  classical  Sommerfeld-like  NRBC.  A  lot  of  work  in  the  mete¬ 
orological  literature  is  based  on  using  H\  with  a  specially  chosen  C\ .  Pearson  [7] 
used  a  special  but  constant  value  of  C\ ,  while  in  the  scheme  devised  by  Orlanski  [8] 
and  in  later  improved  schemes  [9]— [12]  the  C\  changes  dynamically  and  locally  in 
each  time-step  based  on  the  solution  from  the  previous  time-step.  Some  of  the 
limited-area  weather  prediction  codes  used  today  are  based  on  such  schemes,  e.g., 
COAMPS  [13].  See  also  the  recent  papers  [14]— [16]  where  several  such  schemes 
are  compared.  In  a  series  of  papers  [22]— [31] ,  Givoli,  Neta  and  van  Joolen  have 
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demonstrated  the  use  of  high  order  Higdon  NRBC  to  solve  the  shallow  water 
equations  with  advection  and  stratification. 

•  The  condition  Hj  involves  up  to  Jth-order  normal  and  temporal  derivatives.  In 
fact,  it  has  the  form 

didt~^l  =  0  ,  (10) 

3=0 

which  is  obtained  by  expanding  (9). 

•  It  is  easy  to  show  (see  Higdon  [1]  for  a  similar  setting)  that  when  a  wave  of  the 

form  r]  =  impinges  on  the  boundary  T e  where  the  NRBC  Hj  is 

imposed,  the  resulting  reflection  coefficient  is 


*=ri 

3= i 


Cj  -  cx 
Cj  +  cx 


(ii) 


Immediately  we  see  that  if  Cj  =  Cx  for  one  of  the  j’s  then  R  =  0,  namely  there 
is  no  reflection  and  the  NRBC  is  exact.  Moreover,  we  see  that  the  reflection 
coefficient  is  a  product  of  J  factors,  each  of  which  is  smaller  than  1.  This  implies 
that  the  reflection  coefficient  becomes  smaller  as  the  order  J  increases  regardless 
of  the  choice  made  for  the  parameters  Cj.  Of  course,  a  good  choice  for  the  Cj 
would  lead  to  better  accuracy  with  a  lower  order  J,  but  even  if  we  miss  the  correct 
Cj ’s  considerably  (say,  if  we  make  the  simplest  choice  Cj  =  Co  for  j  =  1 , . . . ,  J) , 
we  are  still  guaranteed  to  reduce  the  spurious  reflection  as  we  increase  the  order 
J.  This  is  an  important  property  of  the  Higdon’s  NRBCs  and  is  the  reason  for 
their  robustness. 

•  In  [4] ,  Higdon  points  to  the  possibility  of  a  long-time  instability  that  might  occur 
when  one  uses  a  NRBC  with  high-order  derivatives.  If  the  interior  governing  equa¬ 
tions  and  the  NRBC  both  admit  solutions  at  zero  wave  number  and  frequency, 
and  if  the  data  in  the  problem  include  such  “zero  modes,”  then  a  slowly-growing 
smooth  instability  is  possible.  Whether  this  shows  up  in  practice  depends  on 
the  order  of  the  derivatives  in  the  NRBC  and  the  number  of  spatial  dimensions. 
However,  these  difficulties  do  not  arise  in  the  presence  of  dispersion,  or  if  the  data 
are  confined  to  nontrivial  modes. 


3.  Discretization  of  Higdon’s  NRBCs 


The  Higdon  condition  Hj  is  a  product  of  J  operators  of  the  form  dt  +  Cjdx.  Consider 
the  following  Finite  Difference  (FD)  approximations  (see  e.g.  [17]): 


dt  ~ 


i-sr 

At 


dx 


i-s* 

Ax 


(12) 


In  (12),  At  and  Ax  are,  respectively,  the  time-step  size  and  grid  spacing  in  the  x 
direction,  I  is  the  identity  operator,  and  Sfl  and  S~  are  shift  operators  defined  by 


'-’t  Vpq  Vpq 


n—  1 


c-  r.n  _  n 
'Ipq  'lp—l,q 


(13) 
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Here  and  elsewhere,  is  the  FD  approximation  of  rj(x,y,t)  at  grid  point  (xp,yq)  and 
at  time  tn.  We  use  (12)  in  (9)  to  obtain: 


J 


n 


3= 1 


i-sr 

At 


i-s-\ 

Ax  ) 


vnEq  =  0  • 


(14) 


Here,  the  index  E  corresponds  to  a  grid  point  on  the  boundary  T e-  Higdon  has  solved 
this  difference  equation  (and  also  a  slightly  more  involved  equation  that  is  based  on 
time-  and  space-averaging  approximations  for  dx  and  dt)  for  J  <  3  to  obtain  an  explicit 
formula  for  77^.  This  formula  is  used  to  find  the  current  values  on  the  boundary  T  e 
after  the  solution  in  the  interior  points  and  on  the  other  boundaries  has  been  updated. 
The  formula  for  J  =  2  is  found  in  [6],  and  the  one  for  J  =  3  appears  in  the  appendix 
of  [5].  The  algebraic  complexity  of  these  formulas  increases  rapidly  with  the  order  J. 
It  is  thus  not  surprising  that  we  have  not  found  in  the  literature  any  report  on  the 
implementation  of  the  Higdon  NRBCs  beyond  J  =  3. 

Now  we  show  how  to  implement  the  Higdon  NRBCs  to  any  order  using  a  simple 
algorithm.  To  this  end,  we  first  multiply  (14)  by  At  and  rearrange  to  obtain 


where 


Z  = 


IN 

3= 1 


jl  +  djSt  + 


ejSx 


VEq 


=  0 


CjAt 
a>  ~  1  +  "aT 

dj  =  —1  , 

CjAt 

e  j  = - —  . 

J  Ax 


(15) 


(16) 

(17) 

(18) 


The  coefficient  dj  actually  does  not  depend  on  j.  but  we  keep  this  notation  to  allow 
easy  extensions  to  the  scheme  (see,  e.g.,  [22]).  Now,  this  formula  for  Z  requires  the 
summation  of  3J  terms.  If  we  make  the  simplification 


Cj  =  Cq  V  j  G  1 . . .  J, 


then  our  expansion  for  Z  becomes 


Z 


(aI  +  dSf  +  eSf)Jrffq 

(.55  A*'") 


(19) 


(20) 


where  a  =  J  —  /3  —  7,  and  this  summation  consists  of  only  0+1)0+2)  teraiS)  reducing 
the  computational  time  considerably. 

Note  that  we  need  to  store  77”  values  for  i  =  E,  E  —  1,. . . ,  E  —  J  and  h  =  n,  n  — 

'iq 

1, . . . ,  n  —  J.  In  other  words,  we  have  to  store  the  history  of  the  values  of  y  for  a  layer 
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of  thickness  J  +  1  points  near  the  boundary  T e  and  for  J  +  1  time  levels  (including 
the  current  one).  If  there  are  Ny  grid  points  in  the  y  direction,  then  the  amount  of 
storage  needed  in  a  simple  storage  scheme  is  (J  +  l)2!^.  However,  one  can  save  in 
storage  by  exploiting  the  fact  that  not  all  values  rf}  are  needed,  but  only  those  for 

iq 

which  (E  —  i)  +  (n  —  h)  <  J.  This  is  clear  from  (10)  and  also  from  (15).  For  example, 
the  solution  at  time  tn-j  should  be  stored  only  for  points  on  the  boundary  F  e  itself. 


4.  The  Interior  Scheme 


We  consider  explicit  FD  interior  discretization  schemes  for  the  linearized  Euler  equa¬ 
tions  (1)— (4)  to  be  used  in  conjunction  with  the  Hj  condition.  The  interaction  between 
the  Hj  condition  and  the  interior  scheme  is  a  source  of  concern,  since  simple  choices 
for  an  explicit  interior  scheme  turn  out  to  give  rise  to  long-time  instabilities.  We  have 
tried  the  usual  second-order  centered  difference  scheme 


j]'(x) 


rj(x  +  h)  —  r)(x  —  h) 

2 h 


using  both  Euler’s  method  in  time  and  the  second-order  centered  difference  in  time. 
They  are  stable  for  a  sufficiently  small  time  step  when  used  with  the  boundary  condition 
H i  (which  is  a  Sommerfeld-like  condition  as  previously  mentioned),  but  they  become 
unstable  for  J  >  2.  The  instability  appears  earlier  in  time  when  J  becomes  larger. 
Higdon  [1]  has  proved,  in  the  context  of  the  scalar  Klein-Gordon  equation, 


&tV  ~  C%V2y  +  fr]  =  0  , 


(21) 


that  the  discrete  NRBCs  (14)  are  stable  if  the  interior  scheme  is  the  standard  second- 
order  centered  difference  scheme 


Now  we  shall  show  how  the  linearized  Euler  equations  (1)— (4)  can  be  discretized  in 
such  a  way  as  to  mimic  (22)  and  to  lead  to  a  stable  scheme. 

Let  A t  denote  a  forward  difference  approximation  of  let  V*  denote  a  backward 
difference  approximation  to  the  same,  and  use  similar  notation  for  forward  and  back¬ 
ward  difference  approximations  in  x  and  y.  We  can  then  write  this  wave  equation 
discretization  as 

A tVtulj  =  c2  (A xVxu?j  +  A yVyulj)  (23) 

or 

V tAtufj  =  c2  (yxAxufj  +  V„A yu?j)  (24) 

where  we  use  /  =  0  to  reduce  the  Klein-Gordon  equation  to  the  standard  wave  equation. 
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Try  the  following  discretization  for  (1)— (4): 


V  tP 
A  tu 


A  tv 


V tp 


Po  (V x.U  +  VyV) 
A  xP 


Po 
A  yP 


PO 

~7Po  (Vxu  +  Vyv) 


(25) 


Apply  V x  to  the  second  equation  of  (25),  Vy  to  the  third,  At  to  the  fourth,  and  then 
make  the  appropriate  substitution.  This  gives  us 

A tVtP  =  —  (A XVXP  +  A yS/yp)  ,  (26) 

Po 

exactly  as  desired.  Alternatively,  we  can  switch  the  direction  of  every  spatial  difference 
and  get  the  following  discretization: 


Po  (A xu  +  Ayv) 

Po 

VyP 

Po 

IPo  (A xu  +  A yv) , 

which  is  then  equivalent  to 

A  tVtP  =  —  (' VxAxp  +  VyAyp) 
Po 


V ip  =  - 
A  tu  =  - 

Atv  =  - 
Vf  p  =  - 


(27) 


(28) 


Which  discretization  to  use  is  purely  an  esthetic  decision.  We  have  tried  both,  and  both 
bring  stable,  albeit  asymmetric,  results.  However,  care  must  be  taken  in  programming 
this  scheme  to  ensure  that  the  semi-implicit  computations  of  p  and  p  use  the  correct 
values  for  u  and  v.  In  the  numerical  example  below,  we  use  (25). 


5.  A  Numerical  example 


Let  us  consider  a  simple  numerical  example.  Using  the  discretization  scheme  in  (25),  we 
look  at  a  square  domain  100  km  on  each  side,  subdividing  it  into  a  50  x  50  computational 
domain  with  the  Higdon-like  NRBCs  on  all  four  sides.  Using  a  mean  atmoshperic 
density  of  1.2^|  and  pressure  of  1.01  x  105^  [34],  and  zero  advection,  our  initial 
condition  is  a  biquadratic  pressure  and  density  bulge  in  the  center  of  the  domain: 


P°  = 


Q-21)(30-i)(;;-21)(30-j)\ 
1000  ) 

Po 

(i- 21)(30— i)U- 21)(30— j)\ 

looo  ! 

Po 


21  <  i,j  <  30 
otherwise 

21  <  i,j  <  30 
otherwise 


(29) 
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For  comparison,  our  reference  solution  domain  is  300  km  on  each  side,  with  the  domain 
of  interest  in  the  center.  We  define  the  error  norm  for  each  state  variable  ry  as 


Erj  — 


NrM, 


(30) 


where  Nx,Ny  are  the  number  of  grid  points  in  the  x  and  y  directions,  respectively,  r/j 
is  a  solution  state  variable  using  the  J-order  NRBC,  and  rjo  is  the  reference  solution. 
Our  time  step  is  computed  by 


dt 


\/  dx 2  +  dy 2 

4 Cb 


(31) 


which  is  half  the  CFL  limit,  thus  guaranteeing  stability.  We  run  the  simulation  up 
to  t  =  216,  long  enough  for  the  primary  wave  to  exit  the  computational  domain  with 
the  wave  trough  just  passing  through  the  corners.  The  following  figures  show  the  four 
state  variables  at  the  end  of  the  run  for  J  =  6.  In  each  figure,  the  top  left  shows  the 
computed  solution  using  the  NRBCs,  the  top  right  shows  the  reference  solution,  the 
bottom  left  shows  the  reference  solution  domain  truncated  to  the  size  of  the  computed 
solution’s  domain,  and  the  bottom  right  plots  the  delta  between  the  two  solutions  and 
computes  the  error  norm  as  defined  above. 


Figure  2:  The  solution  for  the  density  p  using  J  =  6 
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X-direction  Velocity 


X-direction  Velocity,  reference  solution 


50  100  150 


X-direction  Velocity,  truncated  reference  solution 


X-direction  Velocity,  absolute  error 


Error  norm  =  0  0  0.00017855 


Figure  3:  The  solution  for  u  using  J  =  6 


Y-direction  Velocity 


Y-direction  Velocity,  reference  solution 


50  100  150 


Figure  4:  The  solution  for  v  using  J 
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J 

p  x  10“5 

MX  10  3 

MX  10  3 

P 

1 

1.7677 

3.7601 

3.7601 

2.0829 

2 

0.43811 

0.94068 

0.94068 

0.51624 

3 

0.18223 

0.38785 

0.38785 

0.21473 

4 

0.10341 

0.22628 

0.22628 

0.12185 

5 

0.083147 

0.16743 

0.16743 

0.097975 

6 

0.068582 

0.17855 

0.17855 

0.080812 

Table  1:  Error  Norms  for  J  G  1,6  with  Discretization  Scheme  (25) 


x  105 
1.02 

1.018 

1.016 

1.014 

1.012 

1.01 

1.008 

1.006 

1.004 

1.002 

50  100  150 

Pressure,  absolute  error 


Pressure,  reference  solution 


Figure  5:  The  solution  for  the  pressure  p  using  J  =  6 

Table  1  shows  the  improvements  as  J  goes  from  1  to  6  (graphs  omitted)  This 
example  demonstrates,  albeit  in  a  simplified  setting,  that  the  linearized  Euler  equations 
are  compatible  with  high-order  Higdon-like  NRBCs. 
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Appendix:  Future  Research 

Here  is  a  list  of  subjects  for  further  investigation  (in  random  order): 
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1.  Thorough  investigation  of  the  numerical  properties  of  the  scheme:  measuring  the 
error  as  a  function  of  the  location  of  the  artificial  boundary;  computing  time  and 
operation  count  as  a  function  of  the  various  parameters  (such  as  J  and  the  number 
of  grid  points  on  the  boundary);  stability  with  various  interior  schemes;  etc. 

2.  Implementing  the  scheme  with  auxiliary  variables,  using  FDs. 

3.  Implementing  the  scheme  with  auxiliary  variables  using  FEs. 

4.  Experimenting  with  the  use  of  the  Higdon  conditions  with  the  Nonlinear  Euler 
equations  in  the  computational  domain.  (Need  to  find  a  stable  interior  scherne- 
NRBC  combination.) 

5.  Applying  the  scheme  in  the  3D  case. 

6.  Extending  the  scheme  to  the  case  of  the  linearized  Euler  equations  with  a  nonzero 
mean  flow  (advection). 
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